#include <time.h>

#include <stdio.h>
#include <math.h>
#include "eno.h"
#include "fluids.h"

#include "fluids.xtn"

extern void init_u();
extern void do_laplacian();
extern void projection_fdiff();
extern void init_parameters();
extern void reset_distance1();
extern void brute_reset_distance();
extern void resolve_triple_point();
extern void calibrate();
extern void uncalibrate();
extern double cvt_density();
extern void get_current_area();
extern void ldiff();
extern void correct_advection();
extern void init_stream_function();
extern void init_marcus();
extern void dump_marcus();
extern void zero_level_set();
extern void stokes_add_wind(); 

#define plotmax (1000)

extern double x1plot[plotmax],y1plot[plotmax],x2plot[plotmax],y2plot[plotmax];
extern int plot_count;

extern double density_below,density_above,density_middle,
              laplace_above,laplace_below,laplace_middle;
extern double power_max;


static double height_min,height_max;
   
double u[imax][jmax],
       v[imax][jmax],
       d1[imax][jmax],
       d2[imax][jmax],
       d3[imax][jmax],
       d4[imax][jmax],d5[imax][jmax],d6[imax][jmax], /* scratch variables */
       d7[imax][jmax],
       yy[imax][jmax],   /* scratch variable in reset1 and invert */
       s[imax][jmax],
       px[imax][jmax],py[imax][jmax],pp[imax][jmax],
       a00[imax][jmax],a01[imax][jmax],a02[imax][jmax],
       a10[imax][jmax],a11[imax][jmax],a12[imax][jmax],
       a20[imax][jmax],a21[imax][jmax],a22[imax][jmax],
       bij[imax][jmax],
       ut[imax][jmax],ut_hld[imax][jmax],u_hld[imax][jmax],
       vt[imax][jmax],vt_hld[imax][jmax],v_hld[imax][jmax],
       d1t[imax][jmax],d1t_hld[imax][jmax],d1_hld[imax][jmax],
       d2t[imax][jmax],d2t_hld[imax][jmax],d2_hld[imax][jmax],
       d3t[imax][jmax],d3t_hld[imax][jmax],d3_hld[imax][jmax],
       st[imax][jmax],st_hld[imax][jmax],s_hld[imax][jmax],
       pcg_result[imax][jmax],
       x[imax],y[jmax],h,delta_t,delta_t_hld;

char file_name[80];
double xmin,ymin,xmax,ymax,
       dump_stepsize,archive_stepsize,
       t_stop,lambda_zero,laplace_coeff,
       boussinesq_coeff,laplace1_coeff,
       density_spread,viscosity_spread,
       surface_tension_level,
       tension_oil_air,tension_oil_water,
       density_water,density_air,
       density_oil,viscosity_water,
       viscosity_air,viscosity_oil,
       current_t,gravity,delta_t_override,
       distance_ignore,farx,fary,farxy,ab_delta_t,
       top_stream_value,bottom_stream_value,rk_time,wind_speed,
       top_stream_value_global,bottom_stream_value_global,
       wind_parameter;

double curvature_indicator[3];

logical use_axisymmetric,project_used,use_symmetric,
        reinit_level,reinit_all,hou_formulation,tension_debug,
        reposition,colortest,gas_burst_under,reinit_test,axipinch,
        use_bmo,free_boundary_flag,brute_force_flag,emad_flag,puckett_flag,
        marcus_flag,freetest_flag,kill_vof_flag,stokes_flag,
        standing_wave_flag,stokes_breaking_flag,stokes_wind_flag,
        stokes_random,exact_flag,curvature_test,stokes_linear_flag,
        terminal_flag,conserve_mass,zalesakcmp,pressure_test,
        use_debug;
int fluid_select,mx,my,mx_type,my_type,feed_back_dnum,order_space,
    ixmin,ixmax,iymin,iymax,diff_type,advance_type,rk_order,
    inner_loop,loop_bound,eno_order,
    mass_type,use_archive;
static int dump_count=0;
int number_time_steps=0;
static double plot_stepsize=0.00;

double stream_error;

/* compute divergence of velocity using central differencing to insure 
   velocity is divergence free 
   input - u,v
   output - value of divergence printed to terminal */
void divergence_check()
{
double total_error,divergence;
int ix,iy;

  total_error=0.0;
  for_all_points_int
    {
     divergence=(gru(ix+1,iy)+gru(ix+1,iy+1)-gru(ix,iy)-gru(ix,iy+1))+
                (grv(ix,iy+1)+grv(ix+1,iy+1)-grv(ix,iy)-grv(ix+1,iy));
     divergence /= (2.0*h);
     total_error += fabs(divergence);
    }
if (use_debug)
  printf("divergence_check%12.7f%12.8f\n",total_error,total_error/(mx*my));
}   /* divergence_check */


void calc_source(maxsource_u,maxsource_v)
double *maxsource_u,*maxsource_v;
{
int ix,iy;
double density,fu,fv;

fu=fv=0.0;

*maxsource_u=*maxsource_v=0.0;

for_domain
 {
  if (solid_body())
   {
   }
  else if (boussinesq_coeff !=0.0)
   {
    density=cvt_density(d1[ix][iy],d2[ix][iy],d3[ix][iy]);
    fv=boussinesq_coeff*density;
   }
  else if (drop_merge())
   {
    /* turn off force when drops reach each other */
    if (current_t<sqrt(2.0/fabs(gravity)))
     if ((d1[ix][iy]<=0.0) || 
         ( (d3[ix][iy]<=0.0)&&(three_phase())&&(feed_back()) )
        )
      fv=(y[iy]>(ymax+ymin)/2.0 ? -fabs(gravity) : fabs(gravity));
   }
  else if ((fluid_select==147)||(fluid_select==47))   /* water splash */
   {
    /* turn off force when drop reaches surface w/vel of 10.0 */
    /* lets try w/slightly higher impact vel 7/3/94 */
    /* lets have drop hit with 0 acceleration 7/3/94 */
    /* lets use same initial velocity as in paper 7/4/94 */
    /* back to normal 7/14/94 */

    if (current_t<0.02)
     fv=-500.0;
    else
     fv=gravity;
   }
  else if ((! drop_jet())&&(! drop_periodic())&&(fluid_select!=70))
   fv=gravity;

  ut[ix][iy]+=fu;
  vt[ix][iy]+=fv;
  *maxsource_u=max(*maxsource_u,fabs(fu));
  *maxsource_v=max(*maxsource_v,fabs(fv));
 }   /* for_domain */
}  /* calc_source */
  
void adjust_time(dnum)
int dnum;
{
int ix,iy;

  for_domain
    {
     if (dnum==1)
      d1[ix][iy]=-min(d2[ix][iy],d3[ix][iy]);
     else if (dnum==2) 
      d2[ix][iy]=-min(d1[ix][iy],d3[ix][iy]);
     else if (dnum==3) 
      d3[ix][iy]=-min(d1[ix][iy],d2[ix][iy]);
    }
}  /* adjust_time */


void matrix_solve1(AA,xx,bb)
double AA[3][3],xx[],bb[];
{
double AA_save[6][6],bb_save[6],
  alpha,hold_value,hold_sum,bb_sum,AA_sum,hold1;
int ix,iy,iz,NN,hold_iy;

  NN=5;
bb_sum=0.0;
for (ix=0;ix<=NN;ix++)
 bb_sum+=fabs(bb[ix]);

AA_sum=0.0;
  for (ix=0;ix<=NN;ix++)
   {
    bb_save[ix]=bb[ix];
    for (iy=0;iy<=NN;iy++)
{
AA_sum+=fabs(AA[ix][iy]);
     AA_save[ix][iy]=AA[ix][iy];
}
   }
  
  for (ix=0;ix<=NN-1;ix++)
   {
    hold_iy=ix;hold_value=AA[ix][ix];
    for (iy=ix+1;iy<=NN;iy++)
     if (fabs(AA[iy][ix])>hold_value)
      {
       hold_iy=iy;
       hold_value=fabs(AA[iy][ix]);
      }

    if (hold_iy!=ix)
     for (iy=ix;iy<=NN;iy++)
      {
       hold_value=AA[ix][iy];
       AA[ix][iy]=AA[hold_iy][iy];
       AA[hold_iy][iy]=hold_value;
      }
    hold_value=bb[ix];
    bb[ix]=bb[hold_iy];
    bb[hold_iy]=hold_value;

    for (iy=ix+1;iy<=NN;iy++)
     {
      alpha=AA[iy][ix]/AA[ix][ix];
      for (iz=ix;iz<=NN;iz++)
       AA[iy][iz]-=alpha*AA[ix][iz];
      bb[iy]-=alpha*bb[ix];
     }
   }   /* gaussian elimination w/partial pivoting */ 

  for (ix=NN;ix>=0;ix--)
   {
    hold_value=bb[ix];
    for (iy=ix+1;iy<=NN;iy++)
     hold_value-=AA[ix][iy]*xx[iy];
    xx[ix]=hold_value/AA[ix][ix];
   }

  hold_value=0.0;
  hold_sum=0.0;
  for (ix=0;ix<=NN;ix++)
   {
    hold_sum+=fabs(xx[ix]);
    hold1=bb_save[ix];
    for (iy=0;iy<=NN;iy++)
     hold1-=AA_save[ix][iy]*xx[iy];
    hold_value+=fabs(hold1);
   }

  printf("matrix solve validation resid,xx sum %e %e\n",
    hold_value,hold_sum);

}    /* matrix_solve1 */

void advance_position()
{
int ix,iy,ix1,iy1,ix2,iy2,ixval,iyval,ip1,
 lenposx[9],lenposy[9];
double xb,yb,xb1,yb1,xb2,yb2,p1,p2,p3,p4,pval,lenval,lensort[9],
  phival[6],xbval[6],ybval[6],AA[6][6],xx[6],bb[6];

 for_all_points
  if (! point_ok)
   {
    d1[ix][iy]=gddd(d1,ix,iy);
    d2t[ix][iy]=guu(d2t,ix,iy);
    d3t[ix][iy]=gvv(d3t,ix,iy);
   }
 for_all_points
  {
   d2[ix][iy]=x[ix];
   d3[ix][iy]=y[iy];
/*
printf("x,y store ix,iy %10.5f%10.5f%5d%5d\n",x[ix],y[iy],ix,iy);
*/
   if (inner_loop==0)
    {
     d1_hld[ix][iy]=d1[ix][iy];
     d2_hld[ix][iy]=d2[ix][iy];
     d3_hld[ix][iy]=d3[ix][iy];

     d2[ix][iy]+=delta_t*d2t[ix][iy];
     d3[ix][iy]+=delta_t*d3t[ix][iy];
     d2t_hld[ix][iy]=d2t[ix][iy];
     d3t_hld[ix][iy]=d3t[ix][iy];
/*
printf("ix,iy,d2t,d3t,d2,d31st %5d%5d%10.5f%10.5f%10.5f%10.5f\n",
 ix,iy,d2t[ix][iy],d3t[ix][iy],d2[ix][iy],d3[ix][iy]);
*/
    }
   else
    {
     d2[ix][iy]+=delta_t*d2t_hld[ix][iy];
if (d2t_hld[ix][iy]!=0.0)
 printf("badix,iy,d2thld %5d%5d%10.5f\n",ix,iy,d2t_hld[ix][iy]);
     d3[ix][iy]+=delta_t*d3t_hld[ix][iy];
if (d3t_hld[ix][iy]!=1.0)
 printf("badix,iy,d3thld %5d%5d%10.5f\n",ix,iy,d3t_hld[ix][iy]);
     xb=d2[ix][iy];yb=d3[ix][iy];
     ixval=(int) ((xb-xmin)/h+0.5);
     iyval=(int) ((yb-ymin)/h+0.5);
     xb1=xmin+(ixval-0.5)*h;xb2=xb1+h;
     yb1=ymin+(iyval-0.5)*h;yb2=yb1+h;

     p1=d2t[ixval][iyval];p2=d2t[ixval+1][iyval];
     p3=d2t[ixval+1][iyval+1];p4=d2t[ixval][iyval+1];
   
     if ((ixval>mx)||(iyval>my))
      {
       p2=guu(d2t,ixval+1,iyval);
       p3=guu(d2t,ixval+1,iyval+1);
       p4=guu(d2t,ixval,iyval+1);
      }
       
     pval=(xb2-xb)*(yb2-yb)*p1+(xb2-xb)*(yb-yb1)*p4+
      (xb-xb1)*(yb2-yb)*p2+(xb-xb1)*(yb-yb1)*p3;
     pval/=sqr(h);
if (pval!=0.0)
 printf("badix,iy,pval %5d%5d%10.5f\n",ix,iy,pval);
     d2[ix][iy]=(d2[ix][iy]+x[ix])/2.0+delta_t*pval/2.0;
/*
printf("ix,iy,xx,d2thld,pval2,d2%5d%5d%10.5f%10.5f%10.5f%10.5f\n",
  ix,iy,x[ix],d2t_hld[ix][iy],pval,d2[ix][iy]);
*/
     p1=d3t[ixval][iyval];p2=d3t[ixval+1][iyval];
     p3=d3t[ixval+1][iyval+1];p4=d3t[ixval][iyval+1];

     if ((ixval>mx)||(iyval>my))
      {
       p2=gvv(d3t,ixval+1,iyval);
       p3=gvv(d3t,ixval+1,iyval+1);
       p4=gvv(d3t,ixval,iyval+1);
      }

     pval=(xb2-xb)*(yb2-yb)*p1+(xb2-xb)*(yb-yb1)*p4+
      (xb-xb1)*(yb2-yb)*p2+(xb-xb1)*(yb-yb1)*p3;
     pval/=sqr(h);
if (pval!=1.0)
 printf("badix,iy,pval1 %5d%5d%10.5f\n",ix,iy,pval);
     d3[ix][iy]=(d3[ix][iy]+y[iy])/2.0+delta_t*pval/2.0;
/*
printf("ix,iy,yy,d3thld,pval3,d3%5d%5d%10.5f%10.5f%10.5f%10.5f\n",
  ix,iy,y[iy],d3t_hld[ix][iy],pval,d3[ix][iy]);
*/
    }
  }   /* initialize our new x,y */

 for_domain
  {
   for (ix1=0;ix1<3;ix1++)
    for (iy1=0;iy1<3;iy1++)
     {
      xb1=d2[ix+ix1-1][iy+iy1-1];
      yb1=d3[ix+ix1-1][iy+iy1-1];
      lenval=sqrt(sqr(xb1-x[ix])+sqr(yb1-y[iy]));
      lensort[3*ix1+iy1]=lenval;
      lenposx[3*ix1+iy1]=ix1;
      lenposy[3*ix1+iy1]=iy1;
     }

   for (ix1=0;ix1<6;ix1++)
    for (iy1=8;iy1>=1+ix1;iy1--)
     if (lensort[iy1]<lensort[iy1-1])
      {
       ip1=lenposx[iy1];
       lenposx[iy1]=lenposx[iy1-1];
       lenposx[iy1-1]=ip1;
       ip1=lenposy[iy1];
       lenposy[iy1]=lenposy[iy1-1];
       lenposy[iy1-1]=ip1;
       p1=lensort[iy1];
       lensort[iy1]=lensort[iy1-1];
       lensort[iy1-1]=p1;
      }
   for (ix1=0;ix1<6;ix1++)
     {
      phival[ix1]=d1_hld[lenposx[ix1]+ix-1][lenposy[ix1]+iy-1];
      xbval[ix1]=d2[lenposx[ix1]+ix-1][lenposy[ix1]+iy-1];
      ybval[ix1]=d3[lenposx[ix1]+ix-1][lenposy[ix1]+iy-1];
     }
   for (ix1=0;ix1<6;ix1++)
     {
      AA[ix1][0]=xbval[ix1];
      AA[ix1][1]=ybval[ix1];
      AA[ix1][2]=1.0;
      AA[ix1][3]=sqr(xbval[ix1]);
      AA[ix1][4]=sqr(ybval[ix1]);
      AA[ix1][5]=xbval[ix1]*ybval[ix1];
      bb[ix1]=phival[ix1];

printf("x,y,phi %10.5f %10.5f %10.5f\n",AA[ix1][0],AA[ix1][1],
 bb[ix1]);

     }
   matrix_solve1(AA,xx,bb);
   d1[ix][iy]=xx[0]*x[ix]+xx[1]*y[iy]+xx[2]+sqr(x[ix])*xx[3]+
    sqr(y[iy])*xx[4]+x[ix]*y[iy]*xx[5]; 
  }    /* looping for_domain */
}   /* advance_position */

void advance_time(w,wt,wt_hld,w_hld)
double w[imax][jmax],wt[imax][jmax],wt_hld[imax][jmax],w_hld[imax][jmax];
{
int ix,iy;
double target;

  for_domain
    {
     if ((number_time_steps==0)&&(advance_type!=2))
      {
       w_hld[ix][iy]=w[ix][iy];
       w[ix][iy]+=delta_t*wt[ix][iy];
       wt_hld[ix][iy]=wt[ix][iy];
      }
     else if (advance_type==0)
      {
       w_hld[ix][iy]=w[ix][iy];
       w[ix][iy]+=delta_t*wt[ix][iy]+
        (sqr(delta_t)/(2.0*delta_t_hld))*(wt[ix][iy]-wt_hld[ix][iy]);
       wt_hld[ix][iy]=wt[ix][iy];
      }
     else if (advance_type==1)
      {
       if (inner_loop==0)
        {
         w_hld[ix][iy]=w[ix][iy];
         w[ix][iy]+=(delta_t/2.0)*wt_hld[ix][iy];
        }
       else
        {
         w[ix][iy]=delta_t*wt[ix][iy]+w_hld[ix][iy];
         wt_hld[ix][iy]=wt[ix][iy];
        }
      }
     else if (advance_type==2)
      {
       if (inner_loop==0)
        {
         w_hld[ix][iy]=w[ix][iy];
         w[ix][iy]+=delta_t*wt[ix][iy];
         wt_hld[ix][iy]=wt[ix][iy];
        }
       else if (rk_order<3)
        {
         w[ix][iy]=(w[ix][iy]+w_hld[ix][iy])/2.0+
          (delta_t/2.0)*wt[ix][iy];
         wt_hld[ix][iy]=wt[ix][iy];
        }
       else if (inner_loop==1)
        {
         w[ix][iy]=(0.25*w[ix][iy]+0.75*w_hld[ix][iy])+
          (0.25*delta_t)*wt[ix][iy];
         wt_hld[ix][iy]=wt[ix][iy];
        }
       else 
        {
         w[ix][iy]=(2.0*w[ix][iy]+1.0*w_hld[ix][iy])/3.0+
          (2.0*delta_t/3.0)*wt[ix][iy];
         wt_hld[ix][iy]=wt[ix][iy];
        }

      }   /* case of advance_type==2 rk */

    }   /* looping all points on grid - point ok */
}    /* advance_time */

#define sum_int(var) \
 ( (r_minus*(-var[0][0]-var[0][2])+r_plus*(-var[2][2]-var[2][0])+ \
    3*(r_minus*var[0][1]+r_plus*var[2][1]+r_val*(var[1][0]+var[1][2]))+ \
    16.0*r_val*var[1][1])*sqr(h)/24.0 )

double compute_avg(ix,iy,piece_var)
int ix,iy;
double piece_var[imax][jmax];
{
int ixc,iyc;
double r_val,r_minus,r_plus,hold[3][3];

  r_val=r_minus=r_plus=1.0;
  if (use_axisymmetric)
   {
    r_val=x[ix];r_minus=r_val-h;r_plus=r_val+h;
   }
  for (ixc=0;ixc<3;ixc++)
   for (iyc=0;iyc<3;iyc++)
    hold[ixc][iyc]=piece_var[ix+ixc-1][iy+iyc-1];
  return sum_int(hold)/sqr(h);
}  /* compute_avg */

void advance_time_dist(w,wt,wt_hld,w_hld)
double w[imax][jmax],wt[imax][jmax],wt_hld[imax][jmax],w_hld[imax][jmax];
{
int ix,iy,it_count;
double avg_val,r_val,error_val;

  advance_time(w,wt,wt_hld,w_hld);
  return;

 /* disable this stuff */

  if ((number_time_steps==0)&&(inner_loop==0))
   {
    for_all_points
    if (! point_ok)
     w[ix][iy]=gddd(w,ix,iy);
    for_domain
     d2_hld[ix][iy]=compute_avg(ix,iy,w);
   }
  advance_time(d2_hld,wt,wt_hld,w_hld);  /* d2_hld is AVERAGED value */
  for_domain
   w[ix][iy]=d2_hld[ix][iy];  /* initial guess */
  for_all_points
   if (! point_ok)
    w[ix][iy]=gddd(w,ix,iy); 

  error_val=0.0; 
  for_domain
   {
    avg_val=compute_avg(ix,iy,w);
    error_val=max(error_val,fabs(d2_hld[ix][iy]-avg_val));
   } 

  it_count=0;
  while ((error_val>tolerance_matrix)&&(it_count<3000))
   {
if (use_debug)
    printf("area error,it %15.9f %5d\n",error_val,it_count);
    for_domain
     {
      r_val=1.0;
      if (use_axisymmetric)
       r_val=x[ix];
      avg_val=compute_avg(ix,iy,w)-(2.0*r_val*w[ix][iy]/3.0);
      w[ix][iy]=(d2_hld[ix][iy]-avg_val)*3.0/(2.0*r_val);
     }

    for_all_points
     if (! point_ok)
      w[ix][iy]=gddd(w,ix,iy);
    error_val=0.0;
    for_domain
     {
      avg_val=compute_avg(ix,iy,w);
      error_val=max(error_val,fabs(d2_hld[ix][iy]-avg_val));
     }
    it_count++;
   }   /* while error too big */
 
}
 
/* reinitialize the level set.
   for 3 immiscible fluids, the three level set functions must 
   be "calibrated" to avoid "void" spots - especially at triple points.
   input d1,d2,d3
   output d1,d2,d3
*/
void modify_density_var()
{
int ix,iy;
double avg_error,test_error,a1,a2;
/*
 if (use_bmo)
  resolve_triple_point();
*/

  if ((! three_phase()) ||
      ((!feed_back())&&(! reinit_all))
     )
   {
    if (brute_force_flag)
     {

      brute_reset_distance(d1,d2,&avg_error,&test_error);
      get_current_area(&a1,d1);
      get_current_area(&a2,d2);
      printf("area d1,d2 error,test error %10.5f%10.5f%10.5f%10.5f\n",
         a1,a2,avg_error,test_error);

      if (avg_error>0.0)
       {
        printf("updating d1 with d2 ....\n");
        for_domain
          d1[ix][iy]=d2[ix][iy];
       }
     }
    else
     reset_distance1(d1,d1_hld,1);
   }  /* single interface above */
  else
   {
    if (feed_back_dnum==-10)
     calibrate();

#define test_reset(dvar,dvar_hld,dnum1) \
  if (feed_back_dnum != dnum1) reset_distance1(dvar,dvar_hld,dnum1)

    test_reset(d1,d1_hld,1);
    test_reset(d2,d2_hld,2);
    test_reset(d3,d3_hld,3);
 
    if (feed_back_dnum==-10)
     uncalibrate();
   }

 if (use_bmo)
  resolve_triple_point();
}

/* get value of stream function at ix,iy.
   compare stream function that is propagated to that which
   is computed using vorticity - they should be the same.
   input pcg_val,s
   output return value
*/
double calc_stream_point(ix,iy)
int ix,iy;
{
double pcg_val;

  if (point1_ok)
   {
    pcg_val=pcg_result[ix][iy];
    stream_error+=fabs(pcg_val-s[ix][iy]);
    return pcg_val;
   }
  else
   return 0.0;
}

/* compute stream function value on the u,v - grid
   (stream function computed on the pressure - grid (offset))
  input ix,iy
  output return value
*/
double calc_stream(ix,iy)
int ix,iy;
{
double stream_avg;

  stream_avg=0.0;
  stream_avg+=calc_stream_point(ix,iy);
  stream_avg+=calc_stream_point(ix-1,iy);
  stream_avg+=calc_stream_point(ix-1,iy-1);
  stream_avg+=calc_stream_point(ix,iy-1);

  return stream_avg/4.0;
}

/* use central differencing to compute vorticity at ix,iy
   input ix,iy
   output return value.
*/
double calc_vorticity(ix,iy)
int ix,iy;
{
double dv_dx,du_dy;

  dv_dx=(gv(ix+1,iy)-gv(ix-1,iy))/h;
  du_dy=(gu(ix,iy+1)-gu(ix,iy-1))/h;
  return (dv_dx-du_dy);
}

/* dump plot data to disk.
   input uu - data
         d_code - data type (e.g. velocity/level set / stream ... )
*/
void dump_data1(uu,d_code)
double uu[imax][jmax];
int d_code;
{
FILE *dump_stream;
char dump_file[25],dump_attr[3];
int ix,iy,low_iy,high_iy;
double out_float,out_area;
int i_close;

  if (d_code==S_DUMP)
   {
    stream_error=0.0;
    init_stream_function();  /* initializes pcg_result */
   }

  if (number_time_steps>0)
   sprintf(dump_attr,"a");
  else
   sprintf(dump_attr,"w");

  sprintf(dump_file,"%s.%03d",file_name,d_code);
  dump_stream=fopen(dump_file,dump_attr);

  fprintf(dump_stream,"%8d\n",mx);

  for_domain
    {
     if (is_distance(d_code))
      out_float=gdt(d_code-D1_DUMP+1,ix,iy);
     else if (d_code==S_DUMP)
      out_float=calc_stream(ix,iy);
     else if (d_code==W_DUMP)
      out_float=calc_vorticity(ix,iy);
     else if (d_code==M_DUMP)
      out_float=sqrt(sqr(u[ix][iy])+sqr(v[ix][iy]));
     else
      out_float=uu[ix][iy];
     fprintf(dump_stream,"%8.4f\n",out_float);
    }

  i_close=fclose(dump_stream);
  if (d_code==S_DUMP)
   printf("************* stream_error %10.5f ***********\n",stream_error);
}   /* dump_data1 */

/* dump plot data */
void dump_data()
{
logical dump_vorticity,dump_stream,dump_velocity,dump_speed,dump_pressure;
int ix,iy;
double out_float;

if (use_debug)
  printf("dumping full grids to file\n");

  dump_vorticity=dump_stream=dump_velocity=dump_speed=dump_pressure=FALSE;

dump_pressure=FALSE;
dump_velocity=TRUE;
if (! use_debug) dump_velocity=FALSE;

#define debug_stream FALSE

  if (debug_stream)
   {
    stream_error=0.0;
    init_stream_function();  /* initializes pcg_result */
    for_domain
     out_float=calc_stream(ix,iy);
    printf("************* stream_error %10.5f ***********\n",stream_error);
   }

  
  dump_data1(d1,D1_DUMP); 
  if (three_phase() )
   {
    dump_data1(d2,D2_DUMP); 
    dump_data1(d3,D3_DUMP); 
   }
  if (dump_velocity)
   {
    dump_data1(u,U_DUMP);
    dump_data1(v,V_DUMP);
   }
  if (dump_pressure)
   dump_data1(pp,P_DUMP);
  if (dump_stream)
   dump_data1(s,S_DUMP);
  if (dump_vorticity)
   dump_data1(u,W_DUMP);
  if (dump_speed)
   dump_data1(u,M_DUMP);
}   /* dump_data */

void dump_pressure(time_count)
int time_count;
{
FILE *pressure_stream;
char pressure_file[85];
int i_close,plot_ctr;
int ix,iy,j1,j2;
double f1,f2,factor;

if (use_debug)
 printf("dumping pressure file \n");

 sprintf(pressure_file,"%s.pressure",file_name);
 if (time_count>0)
  pressure_stream=fopen(pressure_file,"a");
 else
  pressure_stream=fopen(pressure_file,"w");

 if (time_count==0)
  fprintf(pressure_stream,"%5d\n",my+1);
 for (iy=0;iy<=my;iy++)
  fprintf(pressure_stream," %17.10e\n",pp[0][iy]);

 i_close=fclose(pressure_stream);
}   /* dump_pressure */

void dump_stokes(time_count)
int time_count;
{
FILE *stokes_stream;
char stokes_file[85];
int i_close,plot_ctr;
int ix,iy,j1,j2;
double f1,f2,factor;

if (use_debug)
  printf("dumping contour file\n");

  sprintf(stokes_file,"%s.contour",file_name);
  if (time_count>0)
   stokes_stream=fopen(stokes_file,"a");
  else
   stokes_stream=fopen(stokes_file,"w");
 
/* 
  zero_level_set(mx,my,d1);
 
  if (time_count==0)
   fprintf(stokes_stream,"%8.3f %8.3f %8.3f %8.3f\n",
     xmin,ymin,xmax,ymax);

  fprintf(stokes_stream,"%8d\n",plot_count);
  for (plot_ctr=0;plot_ctr<plot_count;plot_ctr++)
   fprintf(stokes_stream,"%8.3f %8.3f %8.3f %8.3f\n",
     x1plot[plot_ctr],y1plot[plot_ctr],x2plot[plot_ctr],y2plot[plot_ctr]);
*/

  for_all_points
   if (! point_ok)
    d1[ix][iy]=gddd(d1,ix,iy);

  for (ix=1;ix<=mx;ix++)
   {
    y1plot[ix]=0.0;
    for (iy=1;iy<=my;iy++)
     {
      j1=iy;
      j2=iy+1;
      if (iy==my)
       j2=iy;
      f1=d1[ix][j1];f2=d1[ix][j2];
      if ( ((f1<0.0)&&(f2>=0.0))||
           ((f1>=0.0)&&(f2<0.0)) )
       y1plot[ix]=(iy-0.5)*h+h*fabs(f1)/(fabs(f1)+fabs(f2));
     }  /* iy */
   }  /* ix */


  if (time_count==0)
   {
    height_min=ymax;
    height_max=ymin;
    fprintf(stokes_stream,"%5d\n",mx);
    for (ix=1;ix<=mx;ix++)
     {
      if (y1plot[ix]<height_min) height_min=y1plot[ix];
      if (y1plot[ix]>height_max) height_max=y1plot[ix];
     }
   }

  factor=2.0*pi/(xmax-xmin);
  if (use_symmetric) factor /= 2.0;
 
  for (ix=1;ix<=mx;ix++)
   fprintf(stokes_stream," %17.10e\n", (y1plot[ix]) );

  i_close=fclose(stokes_stream);
}
 
/* dump current area of bubble to data file file_name.plt */
void dump_area_data()
{
FILE *area_stream,*top_stream,*bot_stream,*mid_stream,*pinch_stream,
     *minor_stream,*major_stream,*pe_stream,*ke_stream,
     *topv_stream,*botv_stream,*midv_stream;
char area_file[85],top_file[85],bot_file[85],mid_file[85],pinch_file[90],
     minor_file[85],major_file[85],pe_file[85],ke_file[85],
     topv_file[85],botv_file[85],midv_file[85];
double out_area,d_ix1,d_ix2,first_iy,second_iy,
       zero_ix,zero_iy,surface_iy,dt,db,dt0,db0,max_rval,current_rval,
       kinetic_energy,potential_energy,r_val,temp_density,
       fraction1,fraction2,vert_vel1,vert_vel2,first_vel,
       second_vel,l1_ke,l2_ke,linf_ke,ke_temp;
int i_close,ix1,ix2,ix,iy,max_ixdiff,max_ix1,max_ix2,
    front_iy,back_iy;
logical first_found,enforce_bubble;

  sprintf(area_file,"%s.plt",file_name);
  sprintf(top_file,"%s.plt.top",file_name);
  sprintf(bot_file,"%s.plt.bot",file_name);
  sprintf(mid_file,"%s.plt.mid",file_name);
  sprintf(topv_file,"%s.plt.topv",file_name);
  sprintf(botv_file,"%s.plt.botv",file_name);
  sprintf(midv_file,"%s.plt.midv",file_name);
  sprintf(pinch_file,"%s.plt.pinch",file_name);
  sprintf(minor_file,"%s.plt.minor",file_name);
  sprintf(major_file,"%s.plt.major",file_name);
  sprintf(pe_file,"%s.plt.pe",file_name);
  sprintf(ke_file,"%s.plt.ke",file_name);

  if (number_time_steps>0)
   {
    top_stream=fopen(top_file,"a");
    bot_stream=fopen(bot_file,"a");
    mid_stream=fopen(mid_file,"a");
    topv_stream=fopen(topv_file,"a");
    botv_stream=fopen(botv_file,"a");
    midv_stream=fopen(midv_file,"a");
    area_stream=fopen(area_file,"a");
    pinch_stream=fopen(pinch_file,"a");
    minor_stream=fopen(minor_file,"a");
    major_stream=fopen(major_file,"a");
    pe_stream=fopen(pe_file,"a");
    ke_stream=fopen(ke_file,"a");
   }
  else
   {
    top_stream=fopen(top_file,"w");
    bot_stream=fopen(bot_file,"w");
    mid_stream=fopen(mid_file,"w");
    topv_stream=fopen(topv_file,"w");
    botv_stream=fopen(botv_file,"w");
    midv_stream=fopen(midv_file,"w");
    area_stream=fopen(area_file,"w");
    pinch_stream=fopen(pinch_file,"w");
    minor_stream=fopen(minor_file,"w");
    major_stream=fopen(major_file,"w");
    pe_stream=fopen(pe_file,"w");
    ke_stream=fopen(ke_file,"w");
   }

  get_current_area(&out_area,d1);

  zero_ix=0.5-xmin/h;
  d_ix1=floor(zero_ix);
  if (fabs(d_ix1-zero_ix)<almost_same)
   d_ix2=d_ix1;
  else 
   d_ix2=d_ix1+1.0;

  ix1=(int) (d_ix1);
  ix2=(int) (d_ix2);
  if (ix1<1) ix1=ix2=1;

  first_found=FALSE;
  first_iy=second_iy=100.0;

  front_iy=my;
  back_iy=0;
  for_domain
    if (gd1(ix,iy)<=0.0)
     {
      if (iy<front_iy)
       front_iy=iy-1;
      if (iy>back_iy)
       back_iy=iy+1;
     }

  max_ixdiff=0;
  max_rval=0.0;

  for_all_y_back
   if (pointy_ok) 
    {
     max_ix2=0;
     for_all_x
      if ((pointx_ok)&&(x[ix]>=0.0))
       {
        if ((gd1(ix,iy)<=0.0)&&(ix>max_ix2))
         max_ix2=ix;
   
        dt=gd1(ix,iy);db=gd1(ix+1,iy);
        if (((dt<0.0)&&(db>=0.0)) ||
            ((dt>=0.0)&&(db<0.0)))
         {
          current_rval=fabs(x[ix]);
          if ((fabs(dt)>almost_same)||(fabs(db)>almost_same))
           current_rval+=h*fabs(dt)/(fabs(dt)+fabs(db));
          max_rval=max(max_rval,current_rval);
         }
       }   /* pointx_ok -- looping horizontally */
     if ((iy>1)&&(iy>front_iy)&&(iy<back_iy))
      if (max_ixdiff<abs(max_ix2-max_ix1))
       max_ixdiff=abs(max_ix2-max_ix1);
     max_ix1=max_ix2;

     dt=gd1(ix2,iy);
     db=gd1(ix2,iy+1);

     if (((dt<0.0) && (db>=0.0)) ||
         ((dt>=0.0) && (db<0.0)))
      {
       dt0=gd1(ix1,iy);
       db0=gd1(ix1,iy+1);
       if (fabs(dt)+fabs(db)<almost_same)
        fraction1=fraction2=0.0;
       else
        {
         fraction1=fabs(dt)/(fabs(dt)+fabs(db)); 
         fraction2=fabs(dt0)/(fabs(dt0)+fabs(db0));
        }
       vert_vel1=gv(ix2,iy)*(1.0-fraction1)+
         gv(ix2,iy+1)*fraction1;
       vert_vel2=gv(ix1,iy)*(1.0-fraction2)+
         gv(ix1,iy+1)*fraction2;
       vert_vel1=(vert_vel1+vert_vel2)/2.0;

       zero_iy=iy+(fraction1+fraction2)/2.0;
       surface_iy=(zero_iy-1)*h+h/2.0+ymin;
       if (! first_found)
        {
         first_iy=surface_iy;
         first_vel=vert_vel1;
        }

       if ((drop_merge()) || (fluid_select==47) || (fluid_select==147))
        enforce_bubble=TRUE;
       else
        enforce_bubble=FALSE;

       if ( (fabs(surface_iy-first_iy)<2.2) || (! enforce_bubble))
        {
         second_iy=surface_iy;
         second_vel=vert_vel1;
        }
       first_found=TRUE;
      }   /* interface was found */
    }    /* looping all vertical points */

  kinetic_energy=potential_energy=0.0;
  l1_ke=l2_ke=linf_ke=0.0;
  ix1=0;
  for_domain
   {
    r_val=1.0;
    if (use_axisymmetric)
     r_val=x[ix];
    temp_density=cvt_density(d1[ix][iy],d2[ix][iy],d3[ix][iy]);
    ke_temp=r_val*temp_density*
      (sqr(u[ix][iy])+sqr(v[ix][iy]));
    kinetic_energy+=ke_temp;
    l1_ke+=ke_temp;
    l2_ke+=sqr(ke_temp);
    linf_ke=max(linf_ke,ke_temp);
    potential_energy+=r_val*y[iy]*fabs(gravity)*temp_density;
    ix1++;
   }
  kinetic_energy /= (ix1);
  potential_energy /= (ix1);
  l1_ke /= (ix1);
  l2_ke=sqrt(l2_ke/ix1);
  printf("ke->l1,l2,linf %12.7f%12.7f%12.7f\n",l1_ke,l2_ke,linf_ke);

if (use_debug)
  printf("dumping area %10.5f\n",out_area);

  fprintf(area_stream,"%12.7f%10.5f\n",current_t,out_area); 
  fprintf(top_stream,"%10.5f%10.5f\n",current_t,first_iy);
  fprintf(bot_stream,"%10.5f%10.5f\n",current_t,second_iy);
  fprintf(mid_stream,"%10.5f%10.5f\n",current_t,(first_iy+second_iy)/2.0);
  fprintf(topv_stream,"%10.5f%10.5f\n",current_t,first_vel);
  fprintf(botv_stream,"%10.5f%10.5f\n",current_t,second_vel);
  fprintf(midv_stream,"%10.5f%10.5f\n",current_t,(first_vel+second_vel)/2.0);
  fprintf(pinch_stream,"%10.5f%10d\n",current_t,max_ixdiff);
  fprintf(minor_stream,"%10.5f%10.5f\n",current_t,fabs(second_iy-first_iy));
  fprintf(major_stream,"%10.5f%10.5f\n",current_t,2*max_rval);
  fprintf(pe_stream,"%10.5f%10.5f\n",current_t,potential_energy);
  fprintf(ke_stream,"%10.5f%10.5f\n",current_t,kinetic_energy);

  i_close=fclose(area_stream);
  i_close=fclose(top_stream);
  i_close=fclose(bot_stream);
  i_close=fclose(mid_stream);
  i_close=fclose(topv_stream);
  i_close=fclose(botv_stream);
  i_close=fclose(midv_stream);
  i_close=fclose(pinch_stream);
  i_close=fclose(minor_stream);
  i_close=fclose(major_stream);
  i_close=fclose(pe_stream);
  i_close=fclose(ke_stream);
}   /* dump_area_data*/

#define archive(xx) fprintf(archive_stream,"%17.10e\n",xx[ix][iy])
#define archive_d(xx) fprintf(archive_stream,"%6d\n",(xx))
#define archive_r(xx) fprintf(archive_stream,"%17.10e\n",(xx))

/* archive data to file file_name.arc */
void archive_data()
{
FILE *archive_stream;
char archive_file[25];
int ix,iy,i_close;

  sprintf(archive_file,"%s.arc",file_name);
  printf("***************** writing to %s *******************\n",archive_file);
  archive_stream = fopen(archive_file,"w");

  for_domain
   {
    archive(u); archive(v);
    archive(d1); archive(d2); archive(d3);
    archive(s);
    archive(px); archive(py); archive(pp);
    archive(ut); archive(ut_hld);
    archive(vt); archive(vt_hld);
    archive(d1t); archive(d1t_hld);
    archive(d2t); archive(d2t_hld);
    archive(d3t); archive(d3t_hld);
    archive(st); archive(st_hld);
    archive(pcg_result);
   }

  archive_r(current_t);
  archive_d(number_time_steps);
  archive_r(delta_t);archive_r(delta_t_hld);archive_r(ab_delta_t);
  archive_r(power_max);
  
  i_close=fclose(archive_stream);
}   /* archive_data */


#define unarchive(xx) fscan_status=fscanf(archive_stream," %le\n",&(xx[ix][iy]))
#define unarchive_d(xx) fscan_status=fscanf(archive_stream," %d\n",&(xx))
#define unarchive_r(xx) fscan_status=fscanf(archive_stream," %le\n",&(xx))

/* read archive data from file_name.arc */
void get_archive()
{
FILE *archive_stream;
char archive_file[25];
int ix,iy,i_close,fscan_status;

  sprintf(archive_file,"%s.arc",file_name);
  printf("***************** reading %s *******************\n",archive_file);
  archive_stream = fopen(archive_file,"r");
  if (archive_stream==NULL)
   printf("cannot find the archive file.\n");

  for_domain
   {
    unarchive(u); unarchive(v);
    unarchive(d1); unarchive(d2); unarchive(d3);
    unarchive(s);
    unarchive(px); unarchive(py); unarchive(pp);
    unarchive(ut); unarchive(ut_hld);
    unarchive(vt); unarchive(vt_hld);
    unarchive(d1t); unarchive(d1t_hld);
    unarchive(d2t); unarchive(d2t_hld);
    unarchive(d3t); unarchive(d3t_hld);
    unarchive(st); unarchive(st_hld);
    unarchive(pcg_result);
   }

  unarchive_r(current_t);
  rk_time=current_t;
  unarchive_d(number_time_steps);
  unarchive_r(delta_t);unarchive_r(delta_t_hld);unarchive_r(ab_delta_t);
  unarchive_r(power_max);
  
  i_close=fclose(archive_stream);
}   /* archive_data */

/* determine whethor time has reached desired increment.
  in_stepsize - frequency to do a specified action
*/
logical time_to_dump(in_stepsize)
double in_stepsize;
{
  if (in_stepsize==0.0)
    return TRUE;
  else if (! (fabs(current_t)<1e6))
    return FALSE;
  else if ( floor(current_t/in_stepsize) !=
       floor((current_t - delta_t)/in_stepsize)
     )
    return TRUE;
  else
    return FALSE;
}

/*
  1.  read from configuration file eno.dat
  2.  read archive data if specified or initialize data using eno_init.c
  3.  advance time
       a. compute convection terms (eno) - conserve.c
       b. compute viscosity (central diff) - elliptic.c
       c. compute surface tension (for computation of cfl cond. only) - tension.c
       d. compute new time step.
       e. project (contains code for doing surface tension) - fdiff.c
       f. re-initialize distance function - reset1.c
       g. advance time with adam bashforth.
*/
void main()
{
int dsize,tbytes,ix,iy;
int ixmaxh,ixminh;
double biggest_cfl,
       tension_u,tension_v,cfl_tension_u,cfl_tension_v,
       test_cfl,cfl_project,
       max_px,max_py,max_force_u,max_force_v,
       max_laplace,max_tension_u,max_tension_v,max_ldiff,
       max_ut,max_vt,
       tension_ratio1,tension_ratio2,tension_ratio3,
       viscous_ratio1,viscous_ratio2,viscous_ratio3,
       tension_cutoff,viscous_cutoff,cfl_ut_cutoff,
       alt_tension_cutoff,
       max_force,force_cutoff,min_time_step,delta_t_start;
int cnt,cfl_ix,cfl_iy;
logical do_lu;
double l1_total,l2_total,local_error,max_local_error;
int ix_loc,iy_loc;
double marcus_stepsize,pressure_stepsize;
clock_t t0,t1,t2,t3;
logical stokes_not_symmetric;

  power_max=0.0;
  marcus_stepsize=0.05;
  pressure_stepsize=0.01;
  current_t=0.0;
  rk_time=0.0;
  number_time_steps=0;
  curvature_test=FALSE;
  project_used = FALSE;
  reinit_level=TRUE;
  reinit_all=TRUE;
  hou_formulation=TRUE;
  tension_debug=FALSE;
  mass_type=2;   /*  1 - pn-pn-1 hn-1  2 - pn-p0  h0   3 - pn-pn-1 h0 */
  conserve_mass=TRUE;
  advance_type=2;     /* 0-adam bashforth 1-alt ab 2- rk */
  rk_order=2;
  eno_order=3;
  diff_type=1;   /* non-conservative (diff_type=2==> stream fn formulation) */

  init_parameters();

/*
  if ((fluid_select==43)||(fluid_select==44))
   mass_type=1;
*/
/* gas burst problems very unstable */

  if (colortest)
   {
    conserve_mass=TRUE;
    mass_type=2;
   }

  stokes_not_symmetric=((stokes_flag)&&
    ((! standing_wave_flag)||(stokes_wind_flag)||(stokes_random)) );

  use_symmetric= TRUE;  /* we check use_axisymmetric flag first */
  if ((periodic_test()==TRUE)||
      ((fluid_select>=80)&&(fluid_select<=83))||
      (pressure_test)||(puckett_flag)||(marcus_flag)||
      (stokes_not_symmetric)|| (dambreak()) || (fluid_select==91) ||
      (rayleigh_taylor()) )  /* 80-zalesak */
   use_symmetric=FALSE;

  dsize=sizeof(double);
  tbytes=dsize*imax*jmax*70;
  printf("allocating tbytes: %10d\n",tbytes);

  for_all_points
   {
    clr(px);clr(py);clr(pp);
    clr(u);clr(v);
    clr(d1);clr(d2);clr(d3);
    clr(s);
    clr(ut);clr(vt);
    clr(d1t);clr(d2t);clr(d3t);
    clr(st);
    clr(ut_hld);clr(vt_hld);
    clr(d1t_hld);clr(d2t_hld);clr(d3t_hld);
    clr(st_hld);
    clr(u_hld);clr(v_hld);
    clr(d1_hld);clr(d2_hld);clr(d3_hld);
    clr(s_hld);
   }
  delta_t_hld=0.0;



  init_u();  /* INITIALIZE u,v,phi HERE */


  kill_vof_flag=TRUE;
  if (diff_type==1)
   kill_vof_flag=FALSE;
   kill_vof_flag=FALSE;

  if (stokes_breaking_flag)
   surface_tension_level*=(xmax-xmin)/(2.0*pi); 
  /* nondim like longuet-higgins and cokelet */

  delta_t_start=delta_t=h/10000.0;  /* first step is simple euler */

  distance_ignore=(density_spread+1008.0)*h;

  printf("symmetric=>%5d reinit_level=>%5d hou_formulation=>%5d\n",
    use_symmetric,reinit_level,hou_formulation);
  printf("tension_debug=>%5d\n",tension_debug);

  printf("gravity=>%10.5f\n",gravity);
  printf("boussinesq=>%10.5f\n",boussinesq_coeff);
  printf("use_bmo=>%5d\n",use_bmo);
  if (brute_force_flag)
   printf("using brute_force reinit- only for 1 interface\n");
  if (! conserve_mass) printf("mass conservation constraint is off \n");
  else if (mass_type==1) printf("fix applied to pn-pn-1 , h'(n-1) \n");
  else if (mass_type==2) printf("fix applied to pn-p0 , h'(0) \n");
  else if (mass_type==3) printf("fix applied to pn-pn-1 , h'(0) \n");
  else printf("invalid mass_type\n");

  if (density_spread != viscosity_spread) 
   printf("WARNING inconsistent density_spread and viscosity_spread\n");
  else
   printf("density_spread and viscosity_spread are the same\n");

  printf("curvature_test %5d\n",curvature_test);
  printf("advance_type %5d\n",advance_type);
  printf("rk_order %5d\n",rk_order);
  printf("eno_order %5d\n",eno_order);
  printf("da,db  ua,ub %12.7f %12.7f %12.7f %12.7f\n",
    density_above,density_below,laplace_above,laplace_below);
  printf("farx,fary,farxy %10.5f %10.5f %10.5f \n",farx,fary,farxy);

  if ((marcus_flag)||(stokes_not_symmetric))
   init_marcus();

  if (use_archive>0)
   {
    project_used=TRUE;
    get_archive();   /* gets d2 back for brute_force case */

    if ((stokes_flag)&&(use_archive==2))
     {
      printf("adding wind, and dumping updated data...\n");
      stokes_add_wind();
      dump_data();
     }
   }
  else
   {
    if (brute_force_flag)
     modify_density_var();   /* initialize d2 */
    else if (colortest)
     for (cnt=0;cnt<1;cnt++)
      modify_density_var();
    printf("writing initial data\n");
    dump_data();
    dump_area_data();
    if ((marcus_flag)||(stokes_not_symmetric))
     dump_marcus(0);
    if (stokes_flag)
     dump_stokes(0);
    if ((fluid_select==47)||(fluid_select==147))
     dump_pressure(0);
   }

#define tension_c(d1,d2,s) ( (s==0.0) ? 1.0e+9 : (d1+d2)/(2.0*s) )

  tension_ratio1=tension_c(density_below,density_above,surface_tension_level);

  if ((three_phase()) && (feed_back()) )
   {
    tension_ratio2=tension_c(density_below,density_middle,tension_oil_air);
    tension_ratio3=tension_c(density_above,density_middle,tension_oil_water);
   }
  else
   tension_ratio2=tension_ratio3=tension_ratio1;

  printf("tension_ratios %10.5f %10.5f %10.5f\n",
   tension_ratio1,tension_ratio2,tension_ratio3);

  tension_ratio1=min(tension_ratio1,tension_ratio2);
  tension_ratio1=min(tension_ratio1,tension_ratio3);

#define viscous_c(d,mu) ( ((mu==0.0) || (laplace_coeff==0.0)) ? \
   1.0e+9 : d/(mu*laplace_coeff) )

  viscous_ratio1=viscous_c(density_below,laplace_below);
  viscous_ratio2=viscous_c(density_above,laplace_above);
  if ((three_phase()) && (feed_back()) )
   viscous_ratio3=viscous_c(density_middle,laplace_middle);
  else
   viscous_ratio3=viscous_ratio2;

  printf("viscous_ratios %10.5f %10.5f %10.5f\n",
   viscous_ratio1,viscous_ratio2,viscous_ratio3);

  viscous_ratio1=min(viscous_ratio1,viscous_ratio2);
  viscous_ratio1=min(viscous_ratio1,viscous_ratio3);

  viscous_cutoff=(3.0/14.0)*sqr(h)*viscous_ratio1;
  tension_cutoff=sqrt(h*tension_ratio1/(4.0*pi))*h;

  printf("viscous_cutoff,tension_cutoff %10.7f %10.7f\n",
    viscous_cutoff,tension_cutoff);
  printf("power_max %10.5f\n",power_max);

/*  **********************************  start of program loop ****** */

  while (current_t<t_stop-0.00001)
   {
    t0=clock();

    loop_bound=2;
    if (advance_type==2)
     loop_bound=rk_order;

    for (inner_loop=0;inner_loop<loop_bound;inner_loop++)
    {

    biggest_cfl=0.0;
    max_px=max_py=max_force_u=max_force_v=0.0;
    max_laplace=max_tension_u=max_tension_v=max_ldiff=0.0;
    cfl_project=0.0;

    for (cnt=0;cnt<3;cnt++)
     curvature_indicator[cnt]=0.0;

    do_lu=FALSE;
    if (number_time_steps==0)
     do_lu=TRUE;
    else if ((inner_loop!=0)||(advance_type!=1))
     do_lu=TRUE;
    if (do_lu)
    {    /* compute L(u) */

    for_all_points
     {
      clr(ut);clr(vt);
      clr(d1t);clr(d2t);clr(d3t);
      clr(st);
     }


    t2=clock();
    if ((!zalesakcmp)&&(!freetest_flag))
     ldiff(&cfl_ix,&cfl_iy,&max_ldiff);
    else
     max_ldiff=0.0;
    t3=clock();
if (use_debug)
    printf("clock ldiff %10ld\n",t3-t2);

    calc_source(&max_force_u,&max_force_v);

    do_laplacian(&max_laplace);
    /* surface tension not added in until the
       projection */

    t2=clock();
if (use_debug)
    printf("clock force+lap. %10ld",t2-t3);

    for_domain
      {
       max_px=abs_max(max_px,px[ix][iy]);
       max_py=abs_max(max_py,py[ix][iy]);
      }   
if (use_debug)
    printf("projection\n");
    projection_fdiff(&cfl_project);
if (use_debug)
    printf("end projection\n");

    if (pressure_test)
     {
      l1_total=l2_total=max_local_error=0.0;
      ix_loc=iy_loc=0;
      for_domain
        {
         local_error=sqrt(sqr(ut[ix][iy])+sqr(vt[ix][iy]));
         l1_total+=local_error;
         l2_total+=sqr(local_error);
         if (local_error>max_local_error)
          {
           ix_loc=ix;iy_loc=iy;
           max_local_error=local_error;
          }
        }  
      l1_total /= (mx*my);
      l2_total /= (mx*my);
      l2_total = sqrt(l2_total);
      printf("linf,l1,l2,mx,my%10.7f%10.7f%10.7f%5d%5d%5d%5d",
        max_local_error,l1_total,l2_total,mx,my,ix_loc,iy_loc);
     }  /* if pressure_test */
    }   /* if ok to compute Lu */

/* *************** begin code to  compute delta_t *****************   */
   
    if (delta_t_override>0.0)
     {
      if (number_time_steps>0)
       delta_t=delta_t_override;
      else if (puckett_flag)
       delta_t=delta_t_start;
      else 
       delta_t=delta_t_override;

      if (current_t+delta_t>t_stop)
       delta_t=t_stop-current_t+0.000000001;
     }
    else
     {
      max_force=sqrt(sqr(max_force_u)+sqr(max_force_v));
      force_cutoff=( (max_force>2.0) ? h/max_force : h/2.0);

      /* assume maximum nondim vel for my problems is 100.0 */
      min_time_step=(h/100.0);
      if ((fluid_select==44)||(fluid_select==43))
       min_time_step=(h/40000.0);

      min_time_step=min(force_cutoff,min_time_step);
      biggest_cfl=( (max_ldiff>2.0) ? h/max_ldiff : h/2.0);
      cfl_ut_cutoff=( (cfl_project>2.0) ? h/cfl_project : h/2.0);

      cfl_ut_cutoff=max(cfl_ut_cutoff,min_time_step);
      biggest_cfl=max(biggest_cfl,min_time_step);

      biggest_cfl=min(biggest_cfl,viscous_cutoff);

      alt_tension_cutoff=tension_cutoff;
/*
      if ( ((fluid_select==44)&&(current_t<0.3))||
           ((fluid_select==43)&&(current_t<0.5))
         )
       alt_tension_cutoff*=(2.0*sqr(h)*sqrt(h));
*/

      biggest_cfl=min(biggest_cfl,alt_tension_cutoff);
      biggest_cfl=min(biggest_cfl,force_cutoff);

      if (inner_loop==0)
       {
        if (number_time_steps==0)
         {
          ab_delta_t=lambda_zero*biggest_cfl;
          delta_t=delta_t_start;
         }
        else
         delta_t=( (advance_type==1) ? ab_delta_t : lambda_zero*biggest_cfl);
  
        if (current_t+delta_t>t_stop)
         delta_t=t_stop-current_t+0.000000001;

        rk_time=current_t+delta_t;
       }
      else if ((inner_loop==1)&&(rk_order==3)&&(advance_type==2))
       rk_time=current_t+delta_t/2.0;
      else
       ab_delta_t=lambda_zero*biggest_cfl;  

     }   /* delta_t_override=0 */
    
if (use_debug)
{
    printf("dx,dt,t,num %10.5f %12.6f %12.6f %5d\n",
      h,delta_t,current_t,number_time_steps);
    printf("min_time_step %12.7f\n",min_time_step);
}
else
 printf("t=%20.10f\n",current_t);
    

/* **************** end code to compute delta_t **************** */

    t2=clock();

    if ((!freetest_flag)&&(! zalesakcmp))
     {
      if (! solid_body())
       {
        advance_time(u,ut,ut_hld,u_hld);
        advance_time(v,vt,vt_hld,v_hld);
        advance_time(s,st,st_hld,s_hld);
       }

      if (do_lu)
       correct_advection(d1,d1t,d1t_hld,d1_hld);

      if (! three_phase())
       {
        if (exact_flag)
         advance_position();  /* only works for runga kutta */
        else
         advance_time_dist(d1,d1t,d1t_hld,d1_hld);
       }
      else
       {
#define advance_d(dval,dvalt,dvalt_hld,dval_hld,dnum1) \
  if (feed_back_dnum!=dnum1) advance_time_dist(dval,dvalt,  \
   dvalt_hld,dval_hld)
      
        advance_d(d1,d1t,d1t_hld,d1_hld,1);
        advance_d(d2,d2t,d2t_hld,d2_hld,2);
        advance_d(d3,d3t,d3t_hld,d3_hld,3);
    
        if (feed_back_dnum>0)
         adjust_time(feed_back_dnum);
       }
      }    /* ! zalesakcmp */
    t3=clock();
if (use_debug)
    printf("time for advance time clock %10ld\n",t3-t2);


    /* adam bashforth or init euler*/
    if (advance_type!=2)  /* no need for past info when using RKF */
     if ((number_time_steps==0)||(advance_type==0))  
      inner_loop++;
    }   /* looping inner_loop */

    /* DO NOT reinitialize after each runga kutta step */
    t2=clock();
    modify_density_var();
    t3=clock();
if (use_debug)
    printf("reinit clock %10ld\n",t3-t2);

    if (! solid_body())
     divergence_check();
    t2=clock();
if (use_debug)
    printf("divergence check clock %10ld\n",t2-t3);

    number_time_steps++;
    current_t += delta_t;
    rk_time=current_t;
    delta_t_hld=delta_t;
    if (time_to_dump(dump_stepsize))
     dump_data();
    if (time_to_dump(archive_stepsize) )
     archive_data();
    if (time_to_dump(plot_stepsize) )
     dump_area_data();
    if (time_to_dump(pressure_stepsize) )
     if ((fluid_select==47)||(fluid_select==147))
      dump_pressure(number_time_steps);

    if (time_to_dump(marcus_stepsize))
     {
      if ((stokes_not_symmetric)||(marcus_flag))
       dump_marcus(number_time_steps);
      if (stokes_flag)
       dump_stokes(number_time_steps);
     }
    t1=clock();
if (use_debug)
{
    printf("output data clock %10ld\n",t1-t2);
    printf("clock time total %10ld\n",t1-t0);
}

   }    /* while not up to the required time */
} /* main */
